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The ropelength of a space curve is usually defined as the quotient of its length by its thickness: the diameter 
of the largest embedded tube around the knot. This idea was extended to space polygons by Eric Rawdon, 
who gave a definition of ropelength in terms of doubly-critical self-distances (local minima or maxima of the 
distance function on pairs of points on the polygon) and a function of the turning angles of the polygon. A naive 
algorithm for finding the doubly-critical self-distances of an n-edge polygon involves comparing each pair of 
edges, and so takes O(n^) time. In this paper, we describe an improved algorithm, based on the notion of 
octrees, which runs in 0{n\ogn) time. The speed of the ropelength computation controls the performance of 
ropelength-minimizing programs such as Rawdon and Piatek's TOROS. An implementation of our algorithm is 
freely available under the GNU Public License. 
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1. INTRODUCTION 

For a curve in 3-space, ropelength is the quotient of 
the length of the curve by its thickness: the diameter of the 
largest embedded tube around the curve. Minimizing rope- 
length is the same as fixing the diameter of the tube and min- 
imizing its length — if the tube is knotted, we are pulling the 
knot tight, and so the minimum ropelength curves in any knot 
type are often called tight knots. Since the problem is such 
a natural one, the definition of thickness has been discovered 
and rediscovered by several authors[l, 3, 14], with the earli- 
est results known (to these authors) on the problem credited to 
Krotenheerdt and Veit in 1976 [12]. 

In the past decade, there has been a great deal of interest in 
exploring the geometry of tight knots; the definition of thick- 
ness has been refined and fully understood [10], it has been 
shown that C^'^ minimizers exist in each knot type[5, 8, 9], 
some minimizing links have been found[5], and a theory of 
ropelength criticality has started to emerge [4, 21]. The de- 
velopment of this theory has been fueled by a steady stream 
of numerical data on ropelength minimizers, from Pieranski's 
original SONO algorithm[15] and Rawdon's TOROS[16], to 
second-generation efforts such as Smutny and Maddocks' 
biarc computations [6, 19] and the RIDGERUNNER project 
of Cantarella, Piatek, and Rawdon. All of these algorithms 
have in their innermost loops a computation of the ropelength 
of a curve in 3-space. 

Intuitively, the thickness of a tube is controlled locally by 
the curvature of the core curve, and globally by the approach 
of "distant" sections of the tube. Rawdon, in his thesis[17], 
defined a radius of curvature for a corner of a polygon. A 
given corner has two circles which are tangent to both incident 
edges and tangent to one of the edges at its center. He proved 
that we can define a sensible polygonal radius of curvature as 
the radius of the smaller of those two circles. 

More precisely: 
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Definition 1. If P„ is a polygonal curve in 'W? with edges 
ei, . . . , e„, and ai is the turning angle of the polygon made 
by edges and ej+i, then let 

minRad(P„)= min 1 I (i) 

^ ' »6i,...,n |2tan(^) ' 2tan(^) J 

where we take e„+i — e\ if the polygon is a closed curve, and 
take i G 1, . . . , 71 — I otherwise. 

Definition 2. Using the distance function on P„ x P„ given 
by D{x, y) — \x ~ y\, we say that a pair xy of Pn (bounding 
the chord xy) is a pair of closest approach q/P„ if it is a non- 
trivial local minimum of the distance function. The length of 
the shortest such chord is denoted POCA(P„) (and we take 
POCA(P„) = oo if no such chord exists). 

Definition 3. We define the thickness o/P„ by 

Thi(P„) = min{2minRad(P„),P0CA(P„)} . (2) 

We note that the value which Rawdon uses in place 
of POCA(P„) in his original definition of polygonal 
thickness[17] is different. In particular, it is always finite. But 
Rawdon reports that the equivalence of the two definitions fol- 
lows from results in an upcoming paper[13]. 

As computing the radius of curvature at a given corner 
only involves the edges incident to that corner, computing 
minRad(P„) requires only 0{n) time. On the other hand, all 
previous efforts to compute thickness have used some variant 
of Algorithm 1 for computing POCA(P„). This algorithm is 
clearly 0{n^). So we have focused our attention on improv- 
ing the POCA(P„) calculation. 

for i — 1 to n do 

for j — i + 1 ton do 

check a and ej for local min chords; 
compare to previous shortest local min chord; 

end 

end 

Algorithm 1: Standard Algorithm for POCA(P„). 

Our algorithm concentrates on reducing the total number of 
edge-edge checks performed by grouping the edges according 
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FIG. 1 : We see the three cases in the proof of Lemma 2. 1 , from left to 
right an edge-edge pair, a vertex-edge pair, and a vertex-vertex pair. 
In the center and right figures we see T" (x) and T'^ (x), and in the 
righthand figure we also see T~ (y). We have not drawn T^{y), but 
it would be colinear with vjVj+i as with T'^{x). 

to their positions in space into a data structure known in com- 
puter graphics as an octree. We will use the octree to optimize 
the inner loop of Algorithm 1, and show that we can isolate 
a constant-size set of candidate Cj 's for any given in time 
O(logn). The new algorithm will then perform 0{n\ogn) 
edge-edge checks, and one octree construction (which will 
also require time 0{n log n)). 

Before continuing, it is reasonable to ask whether such a 
complicated algorithm can be implemented in a way that pro- 
vides a practical advantage over Algorithm 1. We believe 
that our implementation, liboctrope, does. We give per- 
formance data for some test problems in Section 6. And 
more importantly, we invite interested readers to download 
liboctrope and test the code themselves (http : / /ada . 
math . uga . edu / re search /software/oct rope). 

2. EDGE-EDGE CHECKS 

The quantity POCA(P„) is defined to be the smallest non- 
trivial local minimum of the distance function D{x,y) on 
pairs of points on the polygon P„ . To understand it, we first 
make an observation about the nature of these local minima. 

Lemma 2.1. If we orient the curve Pn and let T~ (x), T^{x) 
denote the inward and outward tangent vectors of P„ at x 
( they are different if and only if x is a vertex with nonzero 
turning angle). Every pair xy which locally minimizes D : 

-P« X P„ ^ R has 

T-{x)-{y-x)>Q T+{x)-{y-x)<0 (3) 
T-{y)-{x~y)>Q T+{y) ■ (x ^ y) < 0. (4) 

We note that if x is in the interior of an edge, then the above 
relations force {x) ■ {x — y) — 0. 

Proof. There are three cases: either both x and y are on the 
interior of an edge, one is an edge point and one a vertex, or 
both are vertices, as shown in Figure 1 . At x, the distance from 
y must not decrease to first order as one moves away from 
X in either direction along the curve: a computation verifies 
that this is equivalent to the first line of the statement of the 
Lemma. A similar argument at y completes the proof. □ 

We now make a definition: 
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FIG. 2: The shaded area represents the region of space in which the 
second point y of a locally minimal pair xy can lie when x is on the 
edge e,; or is the vertex Vi. This region consists of the infinite slab 
of parallel planes normal to which pass tlirough e^, together with 
the wedge extending from vertex Vi in the outward direction from the 
vertex. 

Definition 4. The i''' ramp, Ri of a polygonal curve P„ is the 
union of the planes through edge = WiWi+i with normal 
vector Wi+i — Vi, together with the wedge of vectors w defined 
by the inequalities 

{w-v,)■T-{v^)>0 {w-v,)■T+{v^) <Q. (5) 

See Figure 2. 

This leads naturally to the Lemma: 

Lemma 2.2. If xy is a pair of points on Pn which locally 
minimizes D, and x is on the half-open edge ~ {wi+i}, then 
y is in the ramp Ri. 

Proof. If X is in the interior of the edge. Lemma 2.1 implies 
that xy must be perpendicular to e^, and hence that y is in the 
union of normal planes through e^. If a; is at Vi, the inequalities 
above are those of the statement of Lemma 2.1. □ 

It may seem like we have only rephrased Lemma 2.1. In 
fact, we have gained an important geometric insight about the 
problem- for any edge e^, the ramp Ri will probably be very 
close to a thin slab which only intersects the remainder of P„ 
in a few places (see Figure 3). If we can isolate these intersec- 
tions quickly, we can complete the task of finding POCA(P„) 
by a more detailed comparison of these candidates to Ci. 

3. THE OCTREE DATA STRUCTURE 

With the discussion in Section 2, we have reduced the prob- 
lem of identifying edges Cj which may form locally minimal 
pairs with points on edge to the problem of finding which 
edges of P„ intersect e^'s ramp. To do so efficiently, we will 
need a new data structure for P„: the octree[ll]. 

The octree representation of a collection of points in space 
is a tree where each node represents the bounding box of a 
subset of that collection. The eight daughter nodes of a parent 
represent the bounding boxes of subsets of the points in the 
parent box created by dividing that point set in two in each of 
the coordinate directions. The most detailed octree represen- 
tation of a point set has leaf nodes which each contain a single 
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FIG. 3: A typical ramp in a a trefoil of about 90 edges consists of 
a very thin slab which only intersects the remainder of the knot in 
a few places. If we could isolate ramp-knot intersections quickly, it 
would reduce the number of edge-edge checks required to find the 
shortest POCA. 

point but it is common to stop subdividing when the point sets 
are smaller than some fixed number. Figure 4 illustrates three 
levels of this process for a set of points in the plane, while 
Figure 5 shows the resulting tree. 

From the description above, one can observe that it is easy 
to build an octree using the recursive procedure of Algo- 
rithm 2. 

Data : A node of the tree, a corresponding box B, the 

maximum number of points per box m and a list of 
points. 

Result : An octree representation of this point list, 

if the list of points is no longer than m then 

assign these points to this node; 

make this node a leaf of the tree; 

return; 

else 

partition B into 8 child boxes; 
for each child box do 

create a sublist of points intersecting that box; 

recurse if this sublist is nonempty; 

end 

end 

Algorithm 2: One way to Build an Octree 

For an n-point dataset, if one chooses each box partition 
so that no child box contains more than half the total number 
of points in the parent box, the number of levels in this tree 
is less than log2 n, and one expects this algorithm to run in 
O(nlogn) time. However, this algorithm involves a nontriv- 
ial amount of overhead in keeping track of lists of points, and 
making procedure calls. Since we are very concerned with 
the final performance of our implementation, we now present 
a more insightful octree construction algorithm which has the 
same asymptotic time bound of 0{n log n), but is much faster 
than Algorithm 2 in practice. 

To describe the new algorithm, we start with a numbering 
scheme. As we mentioned before, it is conventional to denote 
the upper right quadrant of the plane by the number 1, and 
proceed counterclockwise to the fourth quadrant on the lower 



right. A more natural numbering scheme assigns each quad- 
rant a 2-digit binary number, d^dy, where d^ = has lower 
values of x (the left hand side) and dx — I has higher values 
of X (the right hand side), while dy ^ denotes lower values 
of y (the bottom half), while dy — 1 denotes higher values of 
y (the top half). For octants in 3-space, we could assign three 
digit binary numbers d^dyd^ similarly. 

Now consider the process of quadtree construction again. 
At the first subdivision, we divide the point set in two parts by 
x-coordinate and by y-coordinate. This gives us four groups 
of points, which we can number as above by the 2-digit binary 
numbers d^dy. These groups are the members of the 4 boxes 
in the next level of the tree, as we saw above. 

But there is something else to notice here: If the collection 
of points is sorted by x and by y, the digits d^. and dy for any 
particular point are the most significant binary digits of that 
point's position in the sorted array. Further, if we continue 
to subdivide the points into fourths by x and y, the next pair 
of binary digits associated to each point, d^idy^ will be the 
next pair of binary digits in that point's position in the x and y 
arrays as well. Again, for octrees the situation is similar, but 
we sort by z as well, and create a sequence of 3-digit binary 
(or 1 -digit octal) numbers. 

Continuing this process, we see that each point in the col- 
lection has a unique octal tag generated by interleaving the bi- 
nary digits of its position in the sorted x, y, and z aiTays. This 
tag specifies its position in the octree. Further, if we made a 
least-first traversal of the octree (descending to octants in the 
order of their octal labels), the order in which we would en- 
counter the points would be by increasing octal tags. These 
observations give rise to a new octree-building algorithm: 

Data : A list of points in R^. 

Result : An octree representation of point list. 

Sort the points by x, y, and z coordinates; 

Shuffle binary digits of array positions to create octal tags; 

Sort again by octal tags; 

Build tree from this traversal-ordered list; 

Algorithm 3: A faster octree-building algorithm 

The problem of building a tree from a traversal-ordered list 
of its contents is a standard one in computer science. Our par- 
ticular solution is discussed in some detail in Section 5 below. 
We note that building the octree from the list also has time 
complexity O(nlogn), since every node in the octree must 
be visited, but that this algorithm is still much faster than the 
previous method of octree construction (Algorithm 2). We are 
among many rediscoverers of this method of octree construc- 
tion, which traces its roots to the "linear quadtree" construc- 
tion of Gargantini[7]. 

4. THE CORE OF THE ALGORITHM 

We can now describe our algorithm. Given any e,;, we must 
identify all edges ej which might be part of a shortest POCA 
with Ci. Such edges must obey two conditions: they must 
intersect e^'s ramp, and they must be closer to than the 



FIG. 4: From left to right, these pictures show three stages in the construction of a quadtree representation (the planar version of an octree 
representation) of a set of points. On the left, the bounding box of the entire point set is computed. This is the root of the tree. In the center, we 
see the points divided in two by x and y coordinates, and then grouped by quadrant into four subcoUections, with bounding boxes as shown. 
On the right, we again divide the subcoUections and group into subquadrants. The resulting 3-level tree is shown in Figure 5. 
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FIG. 5: This picture shows the quadtree constructed in Figure 4 in a more familiar form. The boxes are labelled according to the usual 
numbering convention for quadrants of the plane, where the first quadrant is on the top right and numbering proceedings counterclockwise. 
The final numbers show the number of points in each leaf box, and should be compared to the boxes shown in Figure 4 in the right-hand image. 



shortest POCA found so far. Since both conditions can be 
checked for sub-boxes of the octree, we can use them to elim- 
inate groups of e j from consideration before performing edge- 
edge checks. 

In pseudo-code, this is a collection of n calls to the (recur- 
sive) Algorithm 4 (one for each e^). We refer to the entire 
algorithm (minRad computation, octree construction by Al- 
gorithm 3, and calls to Algorithm 4 for each edge) as Octrope. 

Each call to this algorithm might require it to traverse the 
entire depth of the octree before reaching leaf nodes and per- 
forming the edge-edge checks. Yet this depth is bounded 
above by log2 n, so the expected running time for the algo- 
rithm is O(logn). In pathological cases, many or all of the 
boxes may intersect the ramp. If all the boxes intersect all the 
ramps, this algorithm may be asymptotically slower than the 
naive one: we are forced to visit 0{n log n) tree nodes against 
each of n edges, for a total time complexity of 0{n^ logn). 
We have seen Algorithm 1 outperform Algorithm 4 only for 
a particularly bad class of examples: knots formed by con- 
necting vertices chosen at random inside a fixed volume. (See 
Section 6 for details.) 



Data : An octree node, the current minimum POCA length s, 
the maximum number of edges per box m and a ramp 
from Ci. 

Result : All minimum-length POCAs between a and edges in 
this subtree and (perhaps) a smaller value for the current 
minimum POCA length s. 
if this box is within s of Ci then 

if this box intersects the ramp from Ci then 
if this box is a leaf then 

check the (at most m) edges against a; 
if POCAi- of length < s are found then 
update s; 

return list of minimum length POCAs; 

end 

else 

for each nonempty child node do 

I recurse on the child node; 
end 

end 

end 

end 

Algorithm 4: Recursively identifying candidate Cj 's. 



5. IMPLEMENTATION ISSUES 

While being able to replace an 0{'n?) algorithm with one 
which is O(nlogn) will certainly save time /or large enough 
values ofn, there is no guarantee that this will help with prob- 
lems of practical size. Indeed, Algorithm 4 threatens to con- 
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FIG. 6: Our example is a polygonal Hopf-link approximation com- 
posed of two regular pentagons. The lighter lines are the 9 minimum 
length POCAs for which the Octrope algorithm is searching. They 
extend from the midpoints of one side of each interlocked polygon 
to all the midpoints of the edges of the other. 



sume a fair amount of overhead, while Algorithm 1 involves 
only edge-edge checks, which could be coded very efficiently. 
So in this section we turn our attention from the "n log n" to 
its multiplier — from mathematics to program design. In this 
discussion, we'll refer to function names and prototypes from 
our publically available library version of Octrope, which is 
called liboctrope. 

The depth of the octree. Since searching the octree involves 
some overhead, it is to be expected that we will not get the 
best performance from the deepest octree. Rather, we expect 
it to be more efficient to group some number of edges in each 
box and do simple checks between the current edge and the 
edges in an implicated leaf box. 

We implement this by bounding the maximum number of 
levels in the tree by some t and using that bound to calculate 
TO, the maximum number of edges in any leaf box, by the 
formula m — [jT^] (where \r\ is the least integer greater 
than or equal to r). The value of I can be set by the user, 
using the octree_set_levels call or it will default to ^ = 
[| log2 '^1 ' ^ formula at which we arrived empirically. 

A concrete example. We will now trace through our imple- 
mentation of the fast octree construction procedure of Algo- 
rithm 3 for a particular example: a Hopf link where each edge 
is given by a regular pentagon (see Figure 6 and Table 5). 
In this example, ^ would default to 3 and m would then also 
equal 3. 

Sorting edges. The algorithm begins by gathering all of the 
edges into a single n-element array which we call by_oct. 
To avoid double-checking edge pairs later on, we need each 
edge to "belong" to only one of the leaf boxes in our tree. So 
we identify each edge by its midpoint and store that, as well 
as the edge's length and its starting vertex in by_oct. 

As we create by_oct, we also build by_x, by_y and by_z, 
three 7i-element arrays of pointers to the elements of by_oct . 
We then sort these by x, y, and z order The result is shown in 
Table 5. 

We divide by_x, by_z, and by_z into sections of m 



Component 





1 




2 


voo 


= (14.5,20,0) 


VlO 


= (0,0, 14.5) 


Vol 


= (23.5,-7.6,0) 


Vll 


= (0,27.6,23.5) 


V02 


= (0,-24.7,0) 


Vl2 


= (0,44.7, 0) 


V03 


= (-23.5,-7.6,0) 


Vl3 


= (0,27.6, -23.5) 


V04 


= (-14.5,20,0) 


Vi4 


= (0,0,-14.5) 



TABLE I: The approximate vertices of the pentagonal Hopf link 
shown in Figure 6 are given here. We have rounded the numbers 
to the nearest tenth to simplify the table. This does not affect the 
octree-building procedure under discussion, but would change the 
picture shown in the figure above. 



points each (shown in Table 5 by spacing) and walk through 
them, labeling the edges with the binary numbers of the sec- 
tions in which they lie in the following unusual fashion: if 
X = xiX2 ■ ■ ■ xe^u V = Vi-- ■ Vi-i and z = zi • • • are 
the respective box numbers and their binary representations, 
we interleave those bits to produce a single octal number, 
ziyiXiZ2y2 • • ■ ye-ixe-i- This is the octal tag of Section 3 
above ^ We then sort by_oct by that octal tag, as shown in 
Table 5. 

As we discussed in Section 3, the sorted by_oct array is 
in the same order as that of a traversal of the full octree. 

Building the tree. The actual building of the octree can be 
approached in various directions. We could simply use the 
by_oct array with no futher indexing, traversing it with bi- 
nary searches (an approach which saves space at the expense 
of time). On the other hand, if we are to index it, we can build 
our index in a top-down fashion, establishing the root node 
and building out to the leaf nodes. We can build in a bottom- 
up fashion, partitioning off parts of by_oct as leaf nodes and 
collecting them together in groups until we reach a single top 
node. Or we can (and do) use a "sideways" or "limb-by-Iimb" 
approach. We take an array of £ box pointers and on them 
build the "left-hand limb" of the tree, all the way from the 
smallest numbered leaf box down to the root box. Each of the 
boxes knows its first edge in by_oct and how many edges it 
has (which are grouped together thanks to the octal sort). 

Then we walk once through by_oct, watching the octal 
tags. As long as the tag is the same as the one before it, we 
simply increment the count of edges in that box. When it 
changes, we do a binary XOR with the previous tag to see how 
much they differ (that is, which of the octal digits changed). 
That tells how many of the boxes in this "limb" are complete. 
After some cleanup (which may include pruning the "limb") 
we leave those boxes and the create the new ones necessary 
to hold this edge. Figure 7 shows our example tree after this 



To construct octal tags, we take a single pass simultaneously through 
by_x, by_y, and by_z, starting with the second box, which has binary 
tag 00 ■ • ■ 01. As we walk the an'ays, we spread the bits of the box number 
apart (e.g. 1101 1001000001) using a lookup table similar to that of 
Shaffer[18], shift them left 1 bit for y or 2 bits for z, and OR them with the 
tag constructed so far. The tags are thus built up over time and guaranteed 
to be correct only when we reach the end of the pass. 
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bY_x 


by-Y 


bv 




Cos 


(-19,6.2,0) 


C02 ■ 


(-11.75, -16.15,0) 


Ci3 


. \^u, ±o.o, — / 


eo2 


(-11.75, -16.15,0) 


C()i : 


(11.75,-16.15,0) 


Cl2 


: (0,36.15, -11.75) 


eo4 


(0,20,0) 


Ci4 : 


(0,0,0) 


Coo 


: (19,6.2,0) 


eio 


(0, 13.8, 19) 


Cos • 


(-19,6.2,0) 


Coi 




en 


(0,36.15,11.75) 


Coo : 


(19,6.2,0) 


C02 


: (-11.75,-16.15,0) 


ei2 


(0,36.15,-11.75) 


cio : 


(0, 13.8, 19) 


Cos 


: (-19,6.2,0) 


ei3 


(0, 13.8,-19) 


ci3 : 


(0, 13.8, -19) 


C04 


: (0,20,0) 


ei4 


(0,0,0) 


Co4 : 


(0,20,0) 


Cl4 


: (0,0,0) 


eoi 


(11.75, -16.15,0) 


Cii : 


(0,36.15, 11.75) 


Cll 


: (0,36.15, 11.75) 


Coo 


(19,6.2,0) 


ei2 : 


(0,36.15, -11.75) 


ClO 


: (0,13.8, 19) 



TABLE II: In this table, we see the 10 edges of the pentagons in Table 5 sorted by x, y, and z. The edges are sorted by their midpoints, 
and numbered by the index of their first vertices. The spacing reminds us that since m = 3, we are grouping the midpoints by threes when 
constructing boxes. 



edge 


a: -box 


y-hox 


z-hox 


bits 


octal 


decimal 


eo2 








I 


000100 


048 


4 


Cos 





1 


1 


000110 


068 


6 


Coo 


3 


1 





001011 


138 


11 


Coi 


2 





1 


001100 


148 


12 


ei2 


1 


3 





010011 


238 


19 


CIS 


2 


2 





011000 


308 


24 


ClO 


1 


1 


3 


100111 


478 


39 


ei4 


2 





2 


101000 


508 


40 


eo4 





2 


2 


110000 


608 


48 


en 


1 


2 


2 


110001 


6l8 


49 



TABLE III: This table contains the edges with their box numbers in 
the X, y, and z directions, the binary numbers generated by interleav- 
ing the bits of these box numbers, and the corresponding octal tags 
(in octal and decimal). The data is sorted by octal tag, and so appears 
in the same order in which it appears in the by-oct array. In prin- 
ciple, as many as m edges can share an octal tag, which means that 
they occupy the same leaf node of the resulting octree, but this does 
not happen in our example. 

process is complete. 

Searching the tree. We have now created the tree and can 

move into using it. We can check each edge and its ramp 
against the tree, looking for leaf boxes on which to run edge- 
edge checks. Since POCAs are symmetric, we do not ever 
want to compare the same pair of edges twice. To avoid this 
we do the edge-edge check only if the edge in question pre- 
ceeds our chosen edge in by_oct. By so doing, we can elim- 
inate entire boxes because their lowest edges are not in range. 
The improvement gained from technique has been significant. 

6. PERFORMANCE 

We tested our algorithm using a 1.25 Ghz G4 Macintosh 
computer on high-resolution discretizations of trefoil knots 
and on (open) random walks. We compiled our code with 
gcc 3 . 3 and used the -03 option. To make the tests, 
we compared the run times between liboctrope with the 
tree depth set to 1 and with the default tree depth of ^ = 
[| log2 n\ . When the depth is 1, the octree consists of a sin- 



Algorithm 
Standard Octrope Max depth 

Octree levels 1 9 13 

Edge-edge checks 3, 121, 251 32, 033 8189 
Box/ramp checks 51,131 93,187 

Time 1.9 sec 0.16 sec 0.22 sec 

TABLE IV: This table compares the performance of the 
liboctrope library on a 2500-edge random walk at three levels of 
tree depth: 1 (the standard O(n^) algorithm), [| logj n] (Octrope), 
and 13 (the maximum resolution). 



gle box and "'^ edge-edge checks are performed during 
a run. This turns liboctrope into a fairly efficient imple- 
mentation of Algorithm 1 . 

For both of these classes of knots, Octrope was 
much faster than our reference implementation of Algo- 
rithm 1. Figure 8 shows the relative performance of 
the two algorithms on trefoil knots given by 7(0) = 
(( 1 + I cos 3(9) cos 261, (H- 1 cos 36') sin 26*, | sin 36») . In 
fact, Octrope outperformed Algorithm 1 by an even greater 
margin on random walks. 

For trefoils, Algorithm 1 was sometimes faster than 
liboctrope for very small numbers of edges. Figure 9 
shows the performance of both algorithms near the crossover 
point in some detail. 

To understand the effect of varying the number of levels in 
the octree, we also provide data for a 2499-edge random walk 
in Table 6. Here we see a trade-off between edge-edge checks 
and box/ramp checks as the octree resolution increases. In- 
creasing the number of levels in the octree from 9 to 13 cuts 
the number of final edge-edge checks performed by a factor 
of 4, but doubles the number of box/ramp checks. Since the 
box/ramp checks are more computationally expensive, this is 
not a favorable ratio, and the overall execution time increases. 

The data shows that we have been very effective at reducing 
the number of edge-edge checks. On average, Octrope com- 
pares each edge to less than 13 carefully chosen candidates 
when searching for minimum length POCAs. 

It is worth noting that Octrope is not guaranteed to out- 
perform the standard algorithm, even for large numbers of 
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FIG. 7: These two trees show the octree as initially constructed (top) and after our pruning procedure (bottom). This has grouped some edges 
together in single nodes (such as eo2 and eo?,) and deleted some nodes with only one child (such as -823). It is desirable to eliminate extra 
nodes of this kind, since even though we keep the octal tree as compact in memory as possible, each jump from node to node runs the risk of 
straying outside the memory cache of the processor and incurring a delay as more information is loaded from main memory. 
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Number of edges n 

FIG. 8: This plot shows the time (in seconds) required to find the ropelength for trefoils as a function of the number of edges n in the polygonal 
knot. The timings were computed on a 1.25 Ghz Macintosh G4 computer, and represent averages over 10 runs (for times above one second), 
or 100 runs (for times below one second). The data marked with A comes from liboctrope with the default tree depth of ^ = [| logj n] , 
while the data marked with o comes from liboctrope with the tree depth set to 1 to force n(n — 3)/2 edge-edge checks. 
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FIG. 9: This log-log plot shows the time (in seconds) required to find ropelength for trefoils as a function of the number of edges n in the 
polygonal knot near the crossover point where liboctrope becomes faster than Algorithm 1. As in Figure 8, the timings were computed 
on a 1.25 Ghz Macintosh G4 computer, and represent averages over 10 runs (for times above one second), or 100 runs (for times below one 
second). The data marked with a comes from liboctrope with the default tree depth of £ = [| logj n] , while the data marked with o 
comes from liboctrope with the tree depth set to 1 to force n{n — 3)/2 edge-edge checks. The data shows that liboctrope is faster 
than our implementation of Algorithm 1 for trefoil knots with more than about 120 edges. 



edges. For instance, for random knots constructed by choos- 
ing vertices inside a fixed volume, neither ramp-checking nor 
distance checking eliminates a significant number of pairs 
from consideration. The edges are simply too long, and 
pass too close to one another to decide in advance which 
pairs are likely to control thickness. But even in this case, 
liboctrope was only a few times slower than the standard 
algorithm. 

Currently, minimizing the ropelength of 850-edge knots by 
simulated annealing is a relatively taxing task, requiring a few 
weeks of computer time on a standard desktop machine. Our 
timings above show that the ropelength calculations involved 
in that process can be done 5 times faster using liboctrope 
or that calculating the ropelength of a 2000-edge knot with 
liboctrope would take the same time as finding that of 
the 850-edge knot does now. 



7. CONCLUSIONS AND FUTURE DIRECTIONS 

We have given an outline of an improved algorithm for 
computing the ropelength of polygonal space curves in time 
O(nlogn) and contrasted it to the previous standard algo- 
rithm which required time O(n^). We have implemented 
the algorithm efficiently in ANSI C, and given timings which 
show that our algorithm is also much faster in practice than 
previous methods used in the field. 

The increase in speed from using our method should enable 
researchers to consider significantly more complicated knots, 
and to get much higher-resolution data for simpler knots. Both 
of these are valuable goals. It has always been a goal of the 
geometric knot theory community to apply our results to large 
biomolecules such as DNA and proteins. Since these curves 



may involve thousands of vertices, they have been out of the 
reach of tools based on Algorithm 1 . 

However, our methods do not entirely settle the problem of 
fast ropelength computation. As mentioned in the introduc- 
tion, Cantarella, Fu, Kusner, Sullivan, and Wrinkle[4] have 
discovered tiny straight segments in a ropelength-critical sim- 
ple clasp. These segments are a few one-thousands of one unit 
in length out of a total clasp length of about 6 units (a simi- 
lar clasp has been constructed by Starostin[20]). To resolve 
these very small scale phenomena numerically will require 
ropelength-minimized configurations with tens of thousands 
of edges. 

At about 3 seconds per ropelength computation on a stan- 
dard desktop machine, it would simply be untenable to min- 
imize the ropelength of a 20,000-edge knot using our Ubrary 
and simulated annealing on a desktop machine. However, Oc- 
trope parallelizes well, so one could bring supercomputing 
cluster machines to bear on the problem, reducing the time to 
evaluate a configuration to tenths or hundredths of a second. 
This might allow for a long enough cooling schedule to re- 
solve some small-scale phenomena, but there is no guarantee. 

We are hence considering two further approaches to the 
problem: the use of Edelsbrunner's "segment trees" [22] and 
an approach we call the Multiresolution Ropelength Algo- 
rithm. This algorithm is based on the idea that very high res- 
olution knots can be well-approximated by subsampling the 
vertex set. If one keeps track of the distance between the sub- 
sampled knot and the original, one can again eliminate groups 
of edges from edge-edge checking. The potential advantages 
of this scheme are twofold: first, the construction of the cor- 
responding tree is linear in time and second, the subsampled 
knots can be handled by Octrope itself. The disadvantage 
of the multiresolution algorithm is that it will not help with 
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random walks or other very complicated knots, such as large 
protein backbones, as their subsamples will not be close to the 
original curve. 

We would also like to observe that while the discussion 
above is phrased in terms of polygons, the general octree/ramp 
method is equally applicable to other discretization schemes 
for curves, such as biarcs. In that case, the relative speed ad- 
vantage of this algorithm should be greater, since the "edge- 
edge check" for a pair of arcs or spline segments is much 
slower than the edge-edge check for polygonal edges de- 
scribed above. 

In conclusion, we hope that our algorithm and implemen- 
tation will become a standard software component in numer- 
ical investigations of the ropelength problem. If others can 
improve our code, we hope that they will do so, and invite 
them to contact us. We also hope that our public release of 
the library (the first that we know of in geometric knot the- 
ory since Brakke's Evolver[2]) will inspire others in the field 
to contribute from their personal and laboratory collections 
of code to the public domain. Those interested in obtain- 
ing liboctrope can turn to http : //ada .math . uga . 
edu/research/sof tware/octrope/ for further in- 



formation. 



8. ACKNOWLEDGEMENTS 

The authors are grateful to many colleagues, including Her- 
bert Edelsbrunner, Mark Peletier, and John Sullivan, for dis- 
cussions about ropelength and algorithms. The 2002-2003 
VIGRE group in Geometric Knot Theory (in particular Xan- 
der Faber, Chad MuUikin, and Nancy Wrinkle) contributed 
to our understanding of the computational issues surround- 
ing ropelength, and Monica Shaw and Allison Diana, mem- 
bers of the 2003 Summer Undergraduate Research Experi- 
ence, worked on a prototype implementation of the octrope 
algorithm. Michael Piatek and Eric Rawdon served as the 
liboctrope beta-test team, as well as contributing insight 
on efficient code and library design. The authors would 
also like to acknowledge the support of the National Science 
Foundation through the University of Georgia VIGRE grant 
(DMS-00-89927), DMS-99-02397 (to Cantarella), and DMS- 
02-04826 (to Cantarella and Fu). 



[1] Leonard M. Blumenthal and Karl Menger. Studies in geome- 
try. W. H. Freeman and Co., San Francisco, Calif., 1970. (The 
three-point curvature appears on page 320). 

[2] Kenneth A. Brakkc. The surface evolver. Experiment. Math., 
1(2): 141-165, 1992. 

[3] Gregory Buck and Jeremey Orloff. A simple energy function 
for knots. Topology Appl., 6l{3y.205-2U, 1995. 

[4] Jason Cantarella, Joseph H.G. Fu, Robert B. Kusner, John M. 
Sullivan, and Nancy Wrinkle. Criticality for the Gehring link 
problem. arXiv:math.DG/0402212, 2004. 

[5] Jason Cantarella, Robert B. Kusner, and John M. Sullivan. On 
the minimum ropelength of knots and links. Invent. Math., 
150(2):257-286, 2002. 

[6] Mathias Carlen, Ben Laurie, John H. Maddocks, and Jana 
Smutny. Biarcs, global radius of curvature, and the computa- 
tion of ideal knot shapes. In Physical and Numerical Models in 
Knot Theory and Their Application to the Life Sciences. World 
Scientific, 2005. 

[7] Irene Gargantini. An effective way to represent Quadtrees. 
Commun. ACM, 25(12):905-910, 1982. 

[8] O. Gonzalez and R. de la Llave. Existence of ideal knots. J. 
Knot Theory Ramifications, 12(1): 123-133, 2003. 

[9] O. Gonzalez, J. H. Maddocks, F. Schuricht, and H. von der 
Mosel. Global curvature and self-contact of nonlinearly elas- 
tic curves and rods. Calc. Var. Partial Differential Equations, 
14(l):29-68, 2002. 
[10] Oscar Gonzalez and John H. Maddocks. Global curvature, 
thickness, and the ideal shapes of knots. Proa. Natl. Acad. Sci. 
USA, 96(9):4769-4773 (electronic), 1999. 
[11] Chris L. Jackins and Steven L. Tanimoto. Oct-trees and their 
use in representing three-dimensional objects. Camp. Graphics 



and Image Proc, 14:249-270, 1980. 

[12] Otto Krotenheerdt and Sigrid Veit. Zur Theorie massiver 
Knoten. Wiss. Beitn Martin-Luther-Univ. Halle-Wittenberg 
Reihe MMath., 1:61-14, 1976. 

[13] Kenneth C. Millett, Michael Piatek, and Eric Rawdon. Polyg- 
onal knot space near ropelength-minimized knots. Draft copy 
supplied by authors. As of this draft, the result in question is an 
unstated corollary of Theorem 3.6. 

[14] Alexander Nabutovsky. Non-recursive functions, knots "with 
thick ropes", and self-clenching "thick" hyperspheres. Comm. 
Pure Appl. Math., 48(4):38 1-428, 1995. 

[15] Piotr Pieranski. In search of ideal knots. In Ideal knots, vol- 
ume 19 of Sen Knots Everything, pages 20-41. World Sci. Pub- 
lishing, River Edge, NJ, 1998. 

[16] Eric Rawdon. TOROS: Thickness or Ropelength Optimizing 
System. Personal Communication. 

[17] Eric Rawdon. The Thickness of Polygonal Knots. PhD thesis. 
The University of Iowa, 1997. 

[18] Clifford A. Shaffer. Bit interleaving for Quad- or Octrees. In 
Andrew S. Glassner, editor, Graphics Gems, pages 443^47. 
Morgan Kauffman, 1990. 

[19] Jana Smutny and John Maddocks. Approximation of space 
curves with biarcs. 2004. Preprint. 

[20] Eugene L. Starostin. A constructive approach to modelling the 
tight shapes of some linked structures. Proc. Appl. Math. Meek, 
3:479-480, 2003. 

[21] Heiko von der Mosel and Friedemann Schuricht. Characteri- 
zation of ideal knots. Calculus of Variations, 2003. Online. 
DOI: 1 0. 1 007/S00526-003-02 1 6-y 

[22] A. Zomorodian and H. Edelsbrunner. Fast algorithms for box 
intersections. Intemat. J. Comput. Geom. Appl, 12:143-172, 
2002. 



